library(conos)
source('lib.r')
load('Stromal.Nov.RData')
# Rename MSC-3
ano[ano=='MSC-3']='Fibroblasts'
a2=scon$plotGraph(groups=ano,raster=TRUE,plot.na=F,size=0.1,alpha=0.1,font.size = c(5, 5.5),palette=anoS.palf)
a2
ggsave('Stromal.emb.new.v2.pdf',a2,height=4,width=4)
library(cowplot)
features = c('CDH11','NT5E','ENG','CD34','CXCL12','LEPR')#,'LRRC15'
pl=lapply(sn(features), function(x) {print(x);
scon$plotGraph(gene=x, title=x,size=0.1,alpha=0.7,raster = TRUE)})
b=plot_grid(plotlist=pl, nrow = 2)
b
features = c('FAP','FN1','CD44','ITGA5','INHBA','HTRA1')#,'LRRC15'
pl=lapply(sn(features), function(x) {print(x);
scon$plotGraph(gene=x, title=x,size=0.1,alpha=0.7,raster = TRUE)})
b=plot_grid(plotlist=pl, nrow = 2)
b
# 'TNFSF11','TNFRSF11B'
ano = scon$misc$cell.type %>% Toch()
ano = ano[grepl('Peri',ano)] %>% as.factor()
genes=c('KCNJ8','ABCC9','FN1','COL6A1', 'COL6A2', 'COL6A3','AXL','ADIRF','ACTA2','CNN1','MYH11','TAGLN')#,'IL10','CD163','CCL20','CCL4'
exp <- do.call(rbind,lapply(sn(genes),function(gene) conos:::getGeneExpression(scon,gene)))
exp[is.na(exp)]=0
dim(exp)
cn = intersect(names(ano),rownames(scon$embedding))
p=Dotfig(genes,exp[,cn],ano[cn],cols = c("blue","white", "red"))+xlab('')+ylab('')
p
scon$plotGraph(gene='TNFRSF11B')
scon$plotGraph(gene='TNFSF11')
cao$cell.groups = as.factor(ano)
cao$embedding = cao$embedding[names(ano),]
cao$estimateCellDensity()
p0 <- cao$plotEmbedding(color.by='cell.groups', alpha=alpha, size=size,palette=anoM.palf,raster = TRUE,
title='annotation', show.legend=FALSE, font.size=c(3.5,4),plot.na=F)
pl <- cao$plotCellDensity(add.points=TRUE, show.grid=TRUE,contour.conf = "25%" , contours=c('osteoblasts'),
show.cell.groups=FALSE)
plot_grid(p0, plotlist=pl, nrow = 1)
cao$estimateDiffCellDensity(type='wilcox', adjust.pvalues=TRUE, verbose=FALSE, n.permutations=500)
plot_grid(
p0,
cao$plotDiffCellDensity(type='wilcox', title='raw', legend.position=c(0,1),
size=size, alpha=alpha, adjust.pvalues=FALSE,
gradient.range.quantile=0.975, color.range='symmetric'),
cao$plotDiffCellDensity(type='wilcox', title='adj', legend.position=c(0,1),
size=size, alpha=alpha, adjust.pvalues=TRUE),
nrow=1
)
p2= cao$plotDiffCellDensity(type='wilcox', title='Zscore', legend.position=c(0.99,0.14),
size=size, alpha=alpha, adjust.pvalues=FALSE,
gradient.range.quantile=0.975, color.range='symmetric')
p2
b=plot_grid(p0, plotlist=pl,p2, nrow = 2)
b
ggsave('stromal.fig1.pdf',b,height=6.3,width=6)
features=c("RAMP2","PLVAP" ,'ENG',"CD36",'NT5E',"CXCL12" , "LEPR" , "TIMP1" , "COL6A2", "FN1" ,
'DCN','C3','APOD','MFAP4',"MFAP4", "SPP1",'RUNX2','LAPTM5','VAMP8','ACP5',
'RGS5','ACTA2', "THY1",'ADIRF')
library(cowplot)
anoT2 = scon$misc$cell.type %>% as.factor()
cname=names(anoT2)
exp <- do.call(rbind,lapply(sn(features),function(gene) conos:::getGeneExpression(scon,gene)))
exp[is.na(exp)]=0
dim(exp)
cname = intersect(cname,colnames(exp))
exp= exp[rowMeans(exp)>0,]
features=intersect(features,rownames(exp))
p=Dotfig(features,exp[,cname],anoT2[cname],cols = c("blue","white", "red"))+xlab('')+ylab('')
p
ggsave('F4B.marker2.pdf',p,height=7.1,width=5.6)
scon$plotGraph(gene='CSF1')
scon$plotGraph(gene='TGFB1')
scon$plotGraph(gene='CSF1R')
scon$plotGraph(gene='TNFSF12')
scon$plotGraph(gene='TNFRSF11B')
# 'TNFSF11','TNFRSF11B'
ano = scon$misc$cell.type %>% Toch()
ano = ano[grepl('MSC',ano)] %>% as.factor()
genes=c('IL8','TGFB1','VEGFA','VEGFB','COL6A2','COL4A2','COL4A1','COL18A1','COL3A1')#,'IL10','CD163','CCL20','CCL4'
exp <- do.call(rbind,lapply(sn(genes),function(gene) conos:::getGeneExpression(scon,gene)))
exp[is.na(exp)]=0
dim(exp)
cn = intersect(names(ano),rownames(scon$embedding))
#stacked.plot(ano[cn],genes,exp,anoT.pal = NULL,'violin.r.pdf',height = 0.7)
exp= exp[rowMeans(exp)>0,]
features=intersect(genes,rownames(exp))
ano = ano[cn]
ano2 = paste(Toch(ano),Toch(scon$misc$stype)[names(ano)])
names(ano2) = names(ano)
ano2 = as.factor(ano2)
p=Dotfig(genes,exp[,cn],ano2[cn],cols = c("blue","white", "red"))
p
ggsave('dot.genes.pdf',p,height=4.1,width=4.5)
ano2[1:4]
a1= scon$plotGraph(gene='TNFSF11',size=0.2,alpha=0.5,raster = TRUE,title= 'TNFSF11')
a2 = scon$plotGraph(gene='TNFRSF11B',size=0.2,alpha=0.5,raster = TRUE,title= 'TNFRSF11B')
ggsave('TNFSF11.exp.pdf',a1,height=2.2,width=2.1)
ggsave('TNFRSF11B.exp.pdf',a2,height=2.2,width=2.1)
a1= scon$plotGraph(gene='TGFB1',size=0.2,alpha=0.5,raster = TRUE,title= 'TGFB1')
ggsave('TGFB1.exp.pdf',a1,height=2.2,width=2.1)
cn2 = names(ano[ano=='MSC-1']) %>% intersect(.,colnames(exp))
table(scon$misc$stype[cn2])
dim(exp)
modify_vlnplot1('TNFRSF11B',scon$misc$stype[cn2],exp,colp = NULL,height = 0.7,pt.size=0.1)
table(scon$misc$cell.type)
ano = scon$misc$cell.type %>% Toch()
#ano = ano[grepl('MSC',ano)] %>% as.factor()
genes=c('IL8','TGFB1','VEGFA','VEGFB','COL6A2','COL4A2','COL4A1','COL18A1','COL3A1','TNFRSF11A')#,'IL10','CD163','CCL20','CCL4'
genes=c('LRP5', 'ALPL', 'FBN1', 'BGLAP','BMP4')#,'IL10','CD163','CCL20','CCL4'
exp <- do.call(rbind,lapply(sn(genes),function(gene) conos:::getGeneExpression(con0,gene)))
exp[is.na(exp)]=0
dim(exp)
cn2 = names(ano[ano=='osteoblasts']) %>% intersect(.,colnames(exp))
table(scon$misc$stype[cn2])
tmp = scon$misc$stype[cn2] %>% Toch() %>% as.factor()
stacked.plot(tmp,genes,exp,fraction.palette1[levels(tmp)],'osteoblasts.pdf',height = 0.6,pt.size =0.1)
table(scon$misc$cell.type)
ano = scon$misc$cell.type %>% Toch()
#ano = ano[grepl('MSC',ano)] %>% as.factor()
genes=c('IL8','TGFB1','VEGFA','VEGFB','COL6A2','COL4A2','COL4A1','COL18A1','COL3A1','TNFRSF11A')#,'IL10','CD163','CCL20','CCL4'
genes=c('CA2','TCIRG1','CLCN7','OSTM1','PLEKHM1')#,'IL10','CD163','CCL20','CCL4'
exp <- do.call(rbind,lapply(sn(genes),function(gene) conos:::getGeneExpression(con0,gene)))
exp[is.na(exp)]=0
dim(exp)
cn2 = names(ano[ano=='osteoclasts']) %>% intersect(.,colnames(exp))
table(scon$misc$stype[cn2])
tmp = scon$misc$stype[cn2] %>% Toch() %>% as.factor()
stacked.plot(tmp,genes,exp,fraction.palette1[levels(tmp)],'osteoclasts.pdf',height = 0.6,pt.size =0.1)
df=getExp_Sample(exp,'TNFRSF11A',ano[cn2],scon$getDatasetPerCell(),scon$misc$stype[cn2],min.num.cell=5,scale=NULL)
dat2 = data.frame(value = df[['score']], cell = df[['group']], sample = df[['sample']] , Type = df[['fraction']])
library('dplyr')
dat2$Type = as.factor(dat2$Type)
library(ggpubr)
sig=compare_means(value ~ Type, data = dat2,method='t.test') #
sig
siglis=split(sig, seq(nrow(sig)))
pair=lapply(siglis,function(x) as.character(x[,2:3]))
pair
fig = plotMeanMedValuesPerCellType(dat2,type='box', notch = FALSE, show.jitter=TRUE, jitter.alpha=0.2,yline=NULL,plot.theme=theme_classic(),
palette = fraction.palette1,order.x = FALSE,ylab = 'TNFRSF11A expression')+
theme(axis.text.y=element_text(size=10,angle = 60))+
stat_compare_means(comparisons = pair,label = "p.signif",hide.ns=TRUE,size=6,tip.lengt=0.01,bracket.size =0.3,label.y.npc = "bottom")+
ylim(c(min(dat2$value),max(dat2$value)*1.15))
fig
ggsave('TNFRSF11A.osteoclasts.boxplot.pdf',height=2.4,width=1.7)
#drawBoxplot('TNFSF11',dat2,'CTNFSF11 expression',fraction.palette1)
cname = scon$misc$cell.type %>% .[.=='MSC-2'] %>% names()
fig1 = modify_vlnplot1('TGFB1',
scon$misc$stype[cname], exp, colp = NULL, pt.size = 0,stack = FALSE)+theme(legend.position='none')
fig1
cao$estimatePerCellTypeDE(independent.filtering=TRUE, test='DESeq2.Wald')
cao$plotNumberOfDEGenes(pvalue.cutoff=1e-3, p.adjust=FALSE, show.size.dependency=TRUE, font.size=3)
tmp = cao$test.results$de$`MSC-2`
tmp = tmp$res
org <- org.Hs.eg.db
cao$estimateOntology(type="GO", org.db=org, n.top.genes=500)
cao$plotOntologyDistribution(type="GO")
cao$plotOntologySimilarities(type = "GO", genes = "up") # Try also down-regulated
library(clusterProfiler)
cao$plotOntology(type="GO", genes="up", cell.subgroup=x, plot='dot', font.size=8)
names(cao$test.results$de)
res=cao$test.results$de$`osteoclasts`$res
library(EnhancedVolcano)
# create custom key-value pairs for 'high', 'low', 'mid' expression by fold-change
# this can be achieved with nested ifelse statements
keyvals <- ifelse(
res$log2FoldChange <= -2 & res$padj < 0.01, 'royalblue',
ifelse(res$log2FoldChange >= 2 & res$padj < 0.01, 'gold',
'black'))
keyvals[is.na(keyvals)] <- 'black'
names(keyvals)[keyvals == 'gold'] <- 'high'
names(keyvals)[keyvals == 'black'] <- 'mid'
names(keyvals)[keyvals == 'royalblue'] <- 'low'
keyvalssize = keyvals
keyvalssize[keyvalssize == 'gold'] <- 1
keyvalssize[keyvalssize == 'black'] <- 3
keyvalssize[keyvalssize == 'royalblue'] <- 3
keyvalssize = as.numeric(keyvalssize)
keyvalssize <- ifelse(
res$log2FoldChange <= -2 & res$padj < 0.01, 3,
ifelse(res$log2FoldChange >= 2 & res$padj < 0.01, 3,
1))
fig=EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
title=NULL,
subtitle=NULL,
caption=NULL,
selectLab = rownames(res)[which(names(keyvals) %in% c('high', 'low'))],
xlab = bquote(~Log[2]~ 'fold change'),
ylab = bquote(~Log[10]~ 'padj'),
pCutoff = 0.001,
FCcutoff = 2,
pointSize = keyvalssize, #3.5,
labSize = 4.5,
#shape = c(6, 4, 2, 11),
colCustom = keyvals,
colAlpha = 1,
legendPosition = 'none',
legendLabSize = 15,
legendIconSize = 5.0,
drawConnectors = TRUE,
widthConnectors = 1.0,
colConnectors = 'black',
arrowheads = FALSE,
gridlines.major = TRUE,
gridlines.minor = FALSE,
border = 'partial',
borderWidth = 1.5,
borderColour = 'black')
fig
ggsave('Volcano.osteoclasts.pdf',fig,height=4,width=4.7)
save.image(file='Stromal.Nov.RData')
saveRDS(scon$misc$cell.type,'cell.ano.Nov.rds')
saveRDS(scon,'../con.stromal.v2.rds')
scon$plotGraph(gene='TNFSF12')
scon$plotGraph(gene='TNFRSF12A')
scon$plotGraph(gene='VEGFA')
scon$plotGraph(gene='VEGFB')
scon$plotGraph(gene='KDR')
scon$plotGraph(gene='FLT1')
p <- ggplot(na.omit(df1),aes(x=cell,y=pc.of.sample,dodge=Group,fill=Group))+geom_boxplot(notch=FALSE,outlier.shape=NA) + geom_point(position = position_jitterdodge(jitter.width=0.1),color=adjustcolor(1,alpha=0.3),pch=19,size=0.5)+theme_classic()+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_text(angle = 90, hjust = 0.5)) +xlab("") +ylab('Total fraction of stromal cells')+theme(legend.position="top") +
scale_fill_manual(values=fraction.palette1)+ggpubr::stat_compare_means(aes(group = Group),label = "p.signif")
p1 <- ggplot(na.omit(df2),aes(x=cell,y=pc.of.sample,dodge=Group,fill=Group))+geom_boxplot(notch=FALSE,outlier.shape=NA) + geom_point(position = position_jitterdodge(jitter.width=0.1),color=adjustcolor(1,alpha=0.3),pch=19,size=0.5)+theme_classic()+
theme(axis.text.x = element_text(angle = 90, hjust = 1), axis.text.y = element_text(angle = 90, hjust = 0.5)) +xlab("") +ylab('Total fraction of stromal cells')+theme(legend.position="top") +
scale_fill_manual(values=fraction.palette1)+ggpubr::stat_compare_means(aes(group = Group),label = "p.signif")
p
p1
ggsave('fraction1.pdf',p,height=3.7,width=1.7)
ggsave('fraction2.pdf',p1,height=3.7,width=3)
names(con0$samples)
ano = scon$misc$cell.type %>% Toch()
#ano = ano[grepl('MSC',ano)] %>% as.factor()
genes=c('CA2','TCIRG1','CLCN7','OSTM1','PLEKHM1')#,'IL10','CD163','CCL20','CCL4'
genes=c('FB1','ALPL','LRP5')#,'IL10','CD163','CCL20','CCL4'
exp <- do.call(rbind,lapply(sn(genes),function(gene) conos:::getGeneExpression(con0,gene)))
exp[is.na(exp)]=0
dim(exp)
cn2 = names(ano[ano=='MSC-1']) %>% intersect(.,colnames(exp))
table(scon$misc$stype[cn2])
table(ano)
dim(exp)
gene='CA2'
#modify_vlnplot1(gene,scon$misc$stype[cn2],exp,colp = NULL,height = 0.7,pt.size=0.1)
tmp = scon$misc$stype[cn2] %>% Toch() %>% as.factor()
stacked.plot(tmp,genes,exp,fraction.palette1[levels(tmp)],'MSC-1.pdf',height = 0.6,pt.size =0.1)
gene = 'ALPL'
scon$plotGraph(gene=gene)
modify_vlnplot1(gene,scon$misc$stype[cn2],exp,colp = NULL,height = 0.9,pt.size=0.1)